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Within a phonon- assisted resonance level model we develop a self-consistent procedure for calcu- 
lating electron transport currents in molecular junctions with intermediate to strong electron-phonon 
interaction. The scheme takes into account the mutual influence of the electron and phonon subsys- 
tems. It is based on the 2 nd order cumulant expansion, used to express the correlation function of 
the phonon shift generator in terms of the phonon momentum Green function. Equation of motion 
(EOM) method is used to obtain an approximate analog of the Dyson equation for the electron 
and phonon Green functions in the case of many-particle operators present in the Hamiltonian. 
To zero-order it is similar in particular cases (empty or filled bridge level) to approaches proposed 
earlier. The importance of self-consistency in resonance tunneling situations (partially filled bridge 
level) is stressed. We confirm, even for strong vibronic coupling, a previous suggestion concerning 
the absence of phonon sidebands in the current vs. gate voltage plot when the source-drain voltage 
is smalfe 
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I. INTRODUCTION 

Molecular electronics is an active area of research with 
possible technological interest in supplementing currently 
available Si based electronics via further miniaturiza- 
tion of electronic devices^. Experiments on conduction 
in molecular junctions are becoming more common 3-4 . 
Early experiments focused on the absolute conduction 
and on trends such as dependence on wire length, molec- 
ular structure, and temperature. An intriguing issue is 
the role played by molecular nuclear motions. Vibronic 
coupling may lead to rotations, conformational changes, 
atomic rearrangements, and chemical reactions induced 
by the electronic current 6 . It is directly relevant to 
the junction heating problem and can manifest itself in 
polaron-type localization effects and inelastic signals in 
current-voltage spectrai&SiiSiiiiiSii&i^ While weak in- 
elastic signatures are observed in measurements of con- 
ductance with strong electronic coupling between the 
molecule and leads, stronger vibronic peaks with pro- 
nounced phonon sidebands (we use the term phonon 
to describe any vibrational excitation) may be observed 
when this coupling is weak. A full Franck-Condon en- 
velope has been observed with weak molecule-electrode 
interactions at both termini 15 . 

Studies of electron-phonon interaction have a long his- 
tory^, however new points for consideration arise in bi- 
ased current carrying junctions. The interpretation of 
electronic transport in molecular junctions has until re- 
cently has been made largely in the context of multi- 
channel scattering problems 1 '■ 18 > 19 - 2 o ! 2i 4 22 j which disre- 
gards the influence of the contact population (manifesta- 
tions of the Pauli principle blocking of scattering chan- 
nels, change of electronic structure, etc.) on inelastic 
process. Such approaches also disregard the influence 
of the electronic subsystem on the phonon dynamics. A 



systematic framework describing transport phenomena of 
many-particle systems which can take these effects into 
account can be developed based on the non-equilibrium 
Green's function (NEGF) formulation 2 ^ 2 ^. 

Phonon-assisted electron transport in molecular sys- 
tems can be classified by the relative time and energy 
scales of the processes involved. The electron lifetime 
in the junction should be compared to the relevant vi- 
brational frequency^, while the strength of the electron- 
phonon coupling should be judged relative to electronic 
matrix elements (coupling to contacts and/or between 
isolated parts of the molecular system). It is useful 
to consider separately the limits of weak and strong 
electron-phonon coupling. The first corresponds to non- 
resonant phonon-assisted electron tunneling mostly en- 
countered in experiments on inelastic electron tunneling 
spectroscopy (IETS) 7 i 12 i 13 i 14 i 27 . With the development 
and advances in scanning tunneling microscopy (STM) 
and scanning tunneling spectroscopy (STS), IETS has 
proven invaluable as a tool for identifying and charac- 
terizing molecular species within the conduction region. 
The use of Migdal-Eliashberg theory 2 ^* 2 ^ is justified in 
this case, whereupon the lowest non-vanishing (second) 
order perturbation in electron-phonon coupling on the 
Keldysh contour leads to the Born approximation (BA) 
for electron dynamics. This approach using BA or its self- 
consistent (in electron Green function) flavor was used 
in several theoretical studieo 30 i 31 i 32 ' 33 i 34 i 35 i 36 . In a re- 
cent publication^ 7 , we used an advanced version of this 
scheme, self-consistent Born approximation (SCBA) for 
both electron and phonon Green functions, to describe 
features (peaks, dips, line shape, and line width) of the 
IETS signal, d 2 I/d& 2 , as a function of the applied volt- 
age $. Sometimes SCBA is used also in the resonant 
tunneling regim o 35 : 38 . This usage is valid only if electron- 
phonon interaction is weak, so that no essential inelastic 
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features (e.g. polaron formation) can be studied in this 
case. 

The other limit is realized in cases of resonant tunnel- 
ing, which is characterized by longer electron lifetime in 
the junction (though it still may be short relative to the 
characteristic phonon frequency) and stronger effective 
electron-phonon coupling. The perturbative treatment 
breaks down in this case which may result in formation 
of a polaron in the junction. Signatures of resonant tun- 
neling driven by an intermediate molecular ion appear as 
peaks in the first derivative dl /d& and may show phonon 
subbands 39 40 ' 41 . Several theoretical studies of this situ- 
ation in tunneling junctions are available2fl'2i^2^ 3 *^'^. 
Most of them22i2ii22 are based on scattering theory con- 
sideration, others^ 3 ^ 4 *^ are based partially on the NEGF 
methodology. However these works disregard the Fermi 
population in the leads as mentioned above. Another 
approach, the non-equilibrium linked cluster expansion 
(NLCE), is based on generalization of the linked cluster 
expansion to nonequilibrium situations 4 ^. This approach 
takes the contact population into account, but appears to 
be unstable for diagrammatic expansion beyond the low- 
est order. Note also that in all the cases mentioned above 
the phonon subsystem is assumed to remain in thermal 
equilibrium throughout the process. The rate equations 
approach often used in the literature for the case of weak 
coupling to the leads^^i 4 ^ is essentially a quasiclassical 
treatment, having an assumption that the tunneling rate 
is much smaller than decoherence rates on the molecular 
subsystem. While the approach is useful in descri bing 
e.g. C*6o center of mass motion as is done in Ref. I47L 
the intramolecular vibrations we are interested in here 
may have a much longer lifetime, thus making the rate 
equations approach inappropriate. An attempt to gener- 
alize single particle approximation of Refs. EaliUla was 
presented in Ref. 0, where an EOM method was used. 
Another interesting approach uses numerical renormal- 
ization group methodology to study inelastic effects in 
conductance in the linear response regime^. 

In this paper we propose an approximate scheme for 
treating electronic transport in cases involving interme- 
diate to strong electron-phonon coupling in tunnel junc- 
tions. We employ the Keldysh contour based description, 
which treats both electron and phonon degrees of free- 
dom in a self-consistent manner. This approach is close 
in spirit to NLCFj 4 ^ in using a cumulant expansion (when 
finding an approximate expression for phonon shift op- 
erators correlation function in terms of phonon Green 
functions). However unlike NLCE the proposed scheme 
is stable and self-consistent, i.e. the influence of tunnel- 
ing current on the phonon subsystem is taken into ac- 
count. It reduces to the scattering theory results in the 
limit where the molecular bridge energies are far above 
the Fermi energy of the leads and provides a scheme for 
analyzing the effect of the electronic current on the en- 
ergetics in the vibrational space in resonance tunneling 
situations (the issue will be discussed elsewhere) . Deriva- 
tion of equations is based on EOM method, which makes 



our approach similar to Ref. However we go beyond 
it in taking into account renormalization of phonon sub- 
system due to coupling to tunneling electron. Thermal 
relaxation of the molecular phonons, not taken into ac- 
count in 49 , is also introduced in our consideration. 

In the next section we introduce the model and de- 
scribe the approximations made. Section [IIII presents the 
procedure of our self-consistent calculation. In section lTvl 
we report numerical results and compare them to results 
obtained within other approaches. Section IV1 concludes . 



II. MODEL 

We consider a simple resonant-level model with the 
electronic level |0 > coupled to two electrodes left (L) 
and right (R) (each a free electron reservoir at its own 
equilibrium). The electron on the resonant level (elec- 
tronic energy £o) is linearly coupled to a single vibra- 
tional mode (phonon) with frequency loq, henceforth re- 
ferred to as the "primary phonon" . The latter is coupled 
to a phonon bath represented as a set of independent 
harmonic oscillators. The system Hamiltonian is (here 
and below we use h = 1 and e = 1) 

H = Sqc)c+ £kc\c k + (V*4c + h.c.) 

ke{L,R} ke{L,R} 

+ u a d^ d + y^^@P 8 bp + M a Q a Jc + J2UpQ a Qp (1) 

13 

where & (c) are creation (destruction) operators for elec- 
trons on the bridge level, cj, (fife) are corresponding op- 
erators for electronic states in the contacts, a* (a) are 
creation (destruction) operators for the primary phonon, 
and bp (bp) are the corresponding operators for phonon 

states in the thermal (phonon) bath. Q a and Qp are 
phonon displacement operators 

Q a = a + a f Qp = bp + (2) 

The energy parameters M a and Up correspond to the 
vibronic and the vibrational coupling respectively. For 
future reference we also introduce the phonon momentum 
operators 

P a = -i(a-tf) Pp = -i(bp-b^ (3) 

In what follows we will refer to the phonon mode a that is 
directly coupled to the electronic system as the "primary 
phonon" . 

Following previous work on strong electron-phonon in- 
teraction 4 ^^^ we start by applying a small polaron 
(canonical or Lang-Firsov) transformationist to the 
Hamiltonian 

H = e §a He~ §a (4) 
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with 



M, 



OJ 



a /„ + 

- K 



(5) 



Under the additional approximation of neglecting the ef- 
fect of this transformation on the coupling of the primary 
phonon to the thermal phonon bath this leads to 



at; 



H = e &c + ^2 £kC k c k 

ke{L,R} 



k£{L,R} 



E UfiQaQfi 



where 



£o — £o 



A 



Ml, 

LU 



(7) 



A is the electron level shift due to coupling to the primary 
phonon and 



X a = exp 



, M a 

UJQ 



(8) 



is primary phonon shift generator. The Hamiltonian © 
is characterized by the absence of direct electron-phonon 
coupling present in Q. This is replaced by renormal- 
ization of coupling to the contacts. Note that the same 
result can be obtained by repeatedly applying the trans- 
formation analogous to JSJ in the case of weak coupling 
between primary phonon and thermal bath and neglect- 
ing renormalization due to thermal bath phonons (see 
Appendix . 

The Hamiltonian © is our starting point for the cal- 
culation of the steady-state current across t he jun ction, 
using the NEGF expression derived in Refs. 12,^531 



Ik = { J^[Z<(E)G>(E)-Z>(E)G<(E)] 



(9) 



Here Y, K ' are lesser/greater projections of the self- 
energy due to coupling to the contact K (K = L.R) 

£<(£) = if K (E)T K (E) (10) 



if K (E)T K (E) 

f K (E)]T K {E) 



j k\^) — ~n ± - jk\^i\>- k\^i (11) 
with Jk{E) the Fermi distribution in the contact K and 



Y K {E) = 2nJ2 \V k \ 2 5(E-e k ) 



(12) 



keK 



The lesser and greater Green functions in © are 
Fourier transforms to energy space of projections onto 
the real time axis of the electron Green function on the 
Keldysh contour 

G(t 1 ,t 2 ) = -i< T c c(n)ct(r 2 ) > H 

= -i< T c c{Ti)X a (n)c\n)Xl{n) > n (13) 



where the subscripts H and H indicate which Hamilto- 
nian, Q or JHJ) respectively, determines evolution of the 
system, and T c is the contour ordering operator. We use 
the second form and make the usual approximation of 
decoupling electron and phonon dynamics 



G(ti,t 2 ) r* G c (ti,t 2 ) K(t 1 ,t 2 ) 



(14) 



J2 ( v kcicX a + h.c.) where 
(6) 



G c (n,r 2 ) = -i < T c c(n)ct(r 2 ) > s (15) 
K(n,T 2 ) = <T c X a { Tl )Xl{T 2 ) > n (16) 



Below we drop the subscript H keeping in mind that 
Hamiltonian © is the one that determines the time evo- 
lution of the system. Previous studies (see e.g4«*a^) 
stopped here and approximated the functions G c and K 
by the electron Green function in the absence of coupling 
to phonons and by the equilibrium correlation function 
of the phonon shift generator X a , respectively. Here we 
go a step further by 'dressing' these terms in the spirit 
of standard diagrammatic technique 16 . This will yield 
a self-consistent approach for the intermediate to strong 
electron-phonon interaction, the strong coupling analog 
of the self-consistent Born approximation (SCBA) used 
in the weak coupling limi^i. 

We start by expressing the shift generator correlation 
function K in terms of the phonon Green function. One 
can show (see Appendix [B| that in the second order cu- 
mulant expansion this connection is 



K(t 1 ,t 2 ) = exp 



iD 



PaPc 



(r 1 ,r 2 )-<P a 2 > 



} (17) 



where P a is the primary phonon momentum operator de- 
fined in © and 



D 



PaP, 



(n,r 2 ) = -i < T c P a { Tl )P a {T 2 ) > (18) 



is the phonon momentum Green function. < P% >= 
iDp' p (t,t) is time independent in steady-state. 

Next we derive self-consistent equations for the elec- 
tron Green function G c (ti,t 2 ), Ea. H15[) . and the phonon 
momentum Green function Dp a p a (t\, t 2 ), Eq. (|18fl . Since 
we consider situations where the electron-phonon cou- 
pling is strong relative to the coupling between the bridge 
electronic level and the contacts it is reasonable to look 
for an expression of second order in V k . We use an equa- 
tion of motion (EOM) method to obtain expressions for 
these Green functions. Since the Hamiltonian © con- 
tains the exponential operator X, Wick's theorem is not 
applicable so we cannot write the Dyson equations to ex- 
press these Green functions in terms of the corresponding 
self-energies. It is nevertheless possible to obtain the fol- 
lowing approximate coupled integral equations for these 
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Green functions (see Appendices IU1 and lD|l 



electron Green functions in the energy domain 



Dp a p a (r,r')=DP aP (r,T') 



(19) 



+ J d ^J D^ Pa (r, Tl ) U PaPa (ri , r 2 ) D^ Pa (r 2 , r') 
G c {T,r')=G^{Ty) (20) 
+ E / rfr i / rfr 2 G c 0) (r, n) Ec.jf (ti, t 2 ) G<°> (t 2 , t') 

K={L,R} Jc Jc 

where the functions Hp a p a and H C: k are given by 

np a P a (r 1 ,r 2 )=^|^| 2 J Dp^(r 1 ,r 2 )-^ ]T |T4| 2 
/? fce{i,H} 

x g k (T2,Ti)G c (T U T 2 ) < T c X a (n)Xl(T 2 ) > +(n <-> T 2 ) 

(21) 

) = E l^| 2 5fc(n,r 2 ) < T c X (r 2 )lt( Tl ) > 



(22) 



Here K = L,R and is the free electron Green function 
for state k in the contacts. Note that Hp a p a and £ Cj _kt 
play here the same role as self-energies in the Dyson equa- 
tion. Dyson-like equations of this kind are often used as 
approximations when Wick's theorem doesn't work (see 
e.g. ReflU). E qs. CZJ, (O and JUJ can be solved self- 
consistently as described in the next Section. 



D 



(0),r 

Pa Pa 



(E) = 



E - uj + i— 



E + luq + i 



G^' r (E) = \e — e — ^^{E) 



(23) 
(24) 



where we use the wide-band limit (see Ref. |33 for 
discussion) for the phonon retarded self-energy due 
to coupling to thermal bath (the retarded projec- 
tion of the first term on the right in Ea. pil) ') 



Iph 



2TrY,\U \ 2 6(E-u J p) 



(25) 



and where the retarded electron self-energy due to 
coupling to the contacts is taken in the form 



<0),r 



(E) = 



<0),r 
J c : R 

(0) 
K 



2 E- E& + iW-, 



(E) (26) 
(K = L,R) (27) 



.(0) 



K 



with Ej^ and being the center and half width 

of the band respectively and Lj? is the escape 
rate to contact K in the electron wide band limit 
(WjP — > oo relative to all other energy parameters 
of the junction). 



III. SELF-CONSISTENT CALCULATION 
SCHEME 



In what follows we use the term "self-energy" for 
the functions Tlp a p a and ~S c ,k , in analogy to the corre- 
sponding functions that appear in true Dyson equations. 
Eq. (jlTJl . (jT§|l and (|2*U|> provide a self-consistent way to 
solve for the phonon and electron Green functions, D PaPa 
and G c and the corresponding self-energies Iip a p a and 

Since the connection between the correlation function 
of the shift generators and the phonon Green function, 
Ea. (|17|l . is exponential, it is easier to express projections 
of the former function in terms of the latter one in the 
time domain. At the same time, zero-order (no coupling 
between electron and phonon) expressions for the lesser 
and greater projections of Green functions G c and Dp a p a 
are easier to write down in the energy domain. As a re- 
sult, we work in both domains and implement fast Fourier 
transform (FFT) to transform between them. The en- 
suing self-consistent calculation scheme consists of the 
following steps 

1. We start with the zero-order retarded phonon and 



2. The lesser and greater projections of the phonon 
and electron Green functions are obtained using the 
Keldysh equation 25 

D<> a (E) = \D Pa p a (E)\ 2 U<-> a (E) (28) 
Gf>>{E) = \G r c {E)\^<->{E) (29) 

In the first iteration step we use the zero-order 
retarded Green functions (|23() and 124|l with the 
phonon self-energy due to coupling to the thermal 
bath in place of n<' > (£') 



\DPj a (E)fU< p >^ h (E) 



iN(E) 



{E T ^y + { lph /2y 



(30) 



-i[l + N(E)} 



7 P h 



{E±u Q y + { lph /2y 



where N(E) is the Bose distribution of the thermal 
phonon bath. In Ea. H30|) upper (lower) sign corre- 
sponding to lesser (greater) projection. Similarly 
the zero-order (in the vibronic coupling) electron 
self-energy (sum of contributions due to left, L, and 
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right, R, contacts) is used in place of E< !> 

^ < (E)=if K (E)T K (E) 
^<(E) = -i[l-f K (E)]T K (E) 
T K (E) = -2Im[S^( J B)] 



(E) 



(e-e^y + {wPy 



(31) 
(32) 

(33) 



where K — L,R and Jk(E) is the Fermi distribu- 
tion in contact K. The lesser and greater projec- 
tions are then transformed to the time domain. 



3. Utilizing Langreth rules 5 ' 55 and using D p p (t) in 
projections of (jT7Jl one can get lesser and greater 
correlation functions for the shift generator opera- 



tor; 



,56 



<Xl(0)X a {t)> = e^^p.W-^Pa^ )] (34) 
<X a (t)Xl(0)> = e 



,(t=0) 



(35) 



4. Using £<■>(*), Gf>>(t), < Xl(0)X a (t) > and < 
X a (t)Xl(Q) > in projections of the second term 
on the right in Eq. (|21[1 yields the lesser and greater 
phonon self-energies due to coupling to the electron 
in the time domain 



n< el (t) = i\ 2 a <xl(o)x a (t)> 



(36) 



{[^■>(t)]*Gf(t) + Z^<(t) [G>(C} 



n 



PaP, 



<el (t)=i\l<X a (t)Xl(0)> 



(37) 



x{Ei°)'>(t) [Gf(t)Y+[^<(t)]*G>(t)} 



Similarly, projections of Ea. (|22|l lead to lesser and 
greater electron self-energies in the time domain 

£<(*) - ^■<(t)<Xl(0)X a (t)> (38) 
E>(t) = E(°)'>(i)<je o (t)Xt(0)> (39) 

The retarded self-energies in time domain can be 
obtained from the lesser and greater counterparts in 
the usual way E£(i) = 0(t) [£>(*) - £<(t)] In prac- 
tice we use the time domain analog of the Lehmann 
representation^ YT c {t) = e~ St [E>(i) - E<(i)] with 
<5 — > to suppress negative time contributions on 
the FFT grid. The retarded phonon self-energy is 
obtained in the same way. Thus calculated self- 
energies are transformed to the energy domain. 

5. The self-energies obtained in the previous step 
are used to update the Green functions (retarded, 
lesser and greater; phonon and electron) in the en- 
ergy domain. Retarded Green functions are 



D r PaPa (E) = [l/D^J a (E)-n PM (E) 
G r c (E) = [E-eo-KiET 1 



(40) 
(41) 



The lesser and greater Green functions are then 
obtained from the Keldysh equation, Eqs. (|2*8|) and 
(j2HJl. Note that the phonon self-energy there con- 
tains contributions due to both coupling to the 
thermal bath and to the electron, while the elec- 
tron self-energy is a sum of contributions from the 
two contacts dressed by the electron-phonon inter- 
action. The Green functions are then transformed 
to time domain. 

6. The updated Green functions of the previous step 
are used in step |3| closing the self-consistent 
loop. Steps 13151 are repeated until convergence is 
achieved. As a test we use electron population on 
the level 



no = Im [Gf(t = 0)] 



(42) 



Convergence is achieved when the absolute change 
of the level population in two subsequent iterations 
is less than a predefined tolerance. Calculational 
grid is chosen fine enough to support the smallest 
energies and times involved and big enough to span 
over the relevant area. 

The numerical steps described above require repeated 
integrations in time and energy spaces. The numerical 
grids used for these integrations are chosen so as to yield 
converged numerical integrals. Typical grid sizes used 
are of order 0.5 — 2.5 x 10 6 points with energy step of 
order 10~ 4 - 10~ 3 eV. 

After convergence we calculate the current according 
to (JjJJ and lfH|) ■ We also use the converged results to get 
the nonequilibrium electronic density of states 



A(E)= l (G>(E)-G<(E)) 



(43) 



and the nonequilibrium electronic distribution in the 
junction 



f(E)=Im[G<(E)] /A(E) 



IV. RESULTS AND DISCUSSION 



(44) 



Here we present numerical results and compare them to 
those available in the literature. We choose the bands in 
both contacts to be the same, with half width of W [ k ] = 
10 eV positioned in such way that the shifted electronic 
level is in the middle of the band eo = Ejp (K = L,R). 
The band width is taken big enough so that results of the 
calculation do not depend on this choice. The tolerance 
for the self-consistent procedure was 10 -6 . 

We start by presenting results for the equilibrium den- 
sity of states A(E). Figure shows results for relatively 
weak electron-phonon interaction M a ~ T. The parame- 
ters used in this calculation are T = 10 K, = 0.02 eV 
[K = L,R), e = 2 eV, uj = 0.2 eV, M a = 0.063 eV, 
Iph = 0.001 eV. This choice corresponds to reorganiza- 
tion energy of ~ 0.02 eV, Eq.©. Results are presented 
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FIG. 1: Equilibrium DOS for relatively weak electron-phonon 
coupling: self-consistent result (solid line), zero-order result 
(dashed line), uncoupled electron (dotted line). Shown are 
case of filled (a), partially filled (b), and empty (c) electron 
level. See text for parameters. 



for several choices of Fermi energy position: Ep — Eq = 
2 eV (a), eV (b), and —2 eV (c). These cases corre- 
spond to filled, partially filled and empty (hole conduc- 
tion, mixed conduction and electron conduction respec- 
tively) . The solid line represents the result of a full self- 
consistent calculation, the dashed line shows the result 
obtained in the first iteration step (zero-order), and the 
dotted line shows the DOS Aq in the absence of coupling 
to phonons. First, in all the pictures the polaron shift 
of the central (elastic) peak position relative to that of 



A is evident. Second, due to the low temperature and 
the relatively weak electron-phonon coupling only one 
phonon emission peak appears in the plot. This peak is 
symmetric relative to the electron level position (central 
peak) for hole and electron transport (compare Fig. 
and c). In the intermediate regime, Fig. both satel- 
lites appear. In this modest range of electron-phonon 
interaction the zero-order and the self-consistent results 
are very close. Note that the zero-order result for filled 
(Fig. Ok) and empty (Fig. QJ;) levels corresponds to sin- 
gle particle transport. In particular, the dashed line in 
Fig. ^ is exactly the scattering theory result presented 
in Ref. ^2 ( see Fig. 2a there). Note also that, at least in 
the low temperature regime, both the scattering theory 
approach o^ and zero-order approaches within NEGF 
in the way it is implemented in 4 ^* 4 ^ do not describe cor- 
rectly the hole part of single particle transport. 

Phonon absorption and emission sidebands are seen at 
higher temperature and weaker coupling between bridge 
and leads (Figure |2J). The parameters of this calculation 
are T = 300 K, = 0.002 eV (K = L, R), e = 2 eV, 
uj = 0.02 eV, M a = 0.02 eV, j ph = 0.001 eV. The 
choice corresponds to the same reorganization energy of 
~ 0.02 eV. However the electron-phonon coupling here 
is much more pronounced due to weaker coupling to the 
contacts (which implies that the electron-phonon inter- 
action is relatively much stronger in this case). Once 
more Fig. |2K is mirror symmetric of Fig. and a po- 
laronic shift of 0.02 eV between elastic peak in A and 
Aq is observed. The stronger effective electron-phonon 
coupling is manifested in the appearance of five emis- 
sion phonon sidebands and three peaks corresponding to 
absorption and in the fact that the difference between 
zero-order and self-consistent results is much more pro- 
nounced here. Indeed, renormalization due to the in- 
terplay between the electron-phonon interaction and the 
electronic populations in the leads (the Fermi distribu- 
tion in the contacts influences the phonons which in turn 
affect the electron energy distribution on the bridge) may 
change peak heights and positions (see Fig. [2Jt and c) or 
even influence the DOS shape drastically (see Fig. |2b). 
This effect is referred to in Ref. 35 as "floating" of the 
phonon sidebands. The zero-order (dashed) curves in 
Fig. [2b and c can be c omp ared to the corresponding 
curves m Fig. 7 of Ref. |4& Note that our scheme in 
zero-order is similar to the NLCE approach of 46 * (both ap- 
proaches utilize a cumulant expansion). However, while 
the NLCE appears to be problematic when going to the 
second cluster approximation, our scheme is rather sta- 
ble when higher order correlations are included by the 
self-consistent procedure. Moreover, we take the mutual 
influence of the electron and phonon subsystems into ac- 
count rather than assuming phonons in thermal equilib- 
rium, as is done in other work 4 ^ 4 ^ 4 ^ 4 ^. As was already 
mentioned, this self-consistent aspect of the calculation 
may change the results substantially. 

Figure |3| presents the nonequilibrium DOS (dotted 
line) and the distribution function (solid line) obtained 
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E(eV) 

FIG. 2: Equilibrium DOS for relatively strong electron- 
phonon coupling: self-consistent result (solid line) and zero- 
order result (dashed line). Shown are case of filled (a), par- 
tially filled (b), and empty (c) electron level. See text for 
parameters. Dashed vertical line indicates the position of the 
DOS peak in the absence of coupling to phonons. 



from the self-consistent calculation. The distribution 
function in the junction in the absence of coupling to the 
phonon (dashed line) is shown for comparison. Param- 
eters of the calculation are T = 300 K, = 0.02 eV 
(K = L,R), £ = 2 eV, uj = 0.2 eV, M a = 0.2 eV, 
Iph = 0.01 eV (corresponds to reorganization energy of 
~ 0.2 cV). Voltage drop is $ = 1.5 V with E F = 1.8 eV 
and [Ik = Ef±<&/2 (K = L.R). This setup corresponds 
to a mixed (both electron and hole) transport through 




1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 
E(eV) 



FIG. 3: Self-consistent calculation of non-equilibrium DOS 
(dotted line, left vertical axis) and non-equilibrium electron 
energy distribution (solid line, right vertical axis). Shown also 
is the energy distribution for the uncoupled electron (dashed 
line). See text for parameters. 

the junction (analog to graph b in Figs. ^ and Local 
minima in the nonequilibrium population / to the left of 
the elastic peak (local maxima to the right) correspond 
to positions of phonon sidebands in the DOS A. The 
effect is due to outscattering of electrons from the en- 
ergy regions (outscattering of holes from or equivalently 
inscattering of electrons into the energy regions) due to 
phonon emission. One sees that strong electron-phonon 
interaction essentially changes the distribution function. 
Similar results were reported in Ref.^(| (see Fig. 6 there). 

Figure 21 presents self-consistent (solid line) and zero- 
order (dashed line) results for the differential conduc- 
tance dl/dQ as a function of the applied source-drain 
voltage Parameters of the calculation are the same as 
in Fig. |2| except that the calculation is done at low tem- 
perature T = 10 K. Phonon assisted resonant tunneling 
reveals itself in conductance peaks associated with dif- 
ferent vibronic resonances, i.e. electronic levels dressed 
by different numbers of phonon excitations. As was men- 
tioned above, renormalization due to electron-phonon in- 
teraction leads to shift of position and height of phonon 
sideband peaks. The zero-order (dashed) curve can be 
compared to Fig. 6 of Ref. 0. A surprising feature is 
formation of an additional peak at <f> ~ 3.25 V (see in- 
set) in the self-consistent (solid) curve. While the low 
temperature of the calculation makes phonon absorption 
unlikely, we suspect that the effect may be a result of 
phonon absorption due to heating of the phonon subsys- 
tem by electron flux. Clearly such a scenario is possible 
only if coupling to the contacts is treated beyond sec- 
ond order, i.e. within our scheme only the self-consistent 
calculation can yield the effect. 

Figure presents current plotted against the energy 
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(xlO") 




FIG. 4: Differential conductance vs. source-drain voltage. 
Shown are self-consistent (solid line) and zero-order (dashed 
line) results. Enlarged phonon absorption peak in the self- 
consistent result is reproduced in the inset. See text for de- 
tails. 



CxHT) 




-0.1 0.1 

e -A-E F (eV) 



FIG. 5: Current vs. gate potential when source-drain volt- 
age is fixed at $ = 0.01 V. The position of the shifted level 
is given relative to the unbiased Fermi energy. Shown are 
the self-consistent (solid line), zero-order (dashed line), and 
uncoupled electron (dotted line) results. See text for details. 



of the bridge electronic level (that can be controlled by 
a gate voltage) at small fixed source-drain voltage. The 
horizontal axis represents the position of the shifted elec- 
tron level relative to the original (unbiased) Fermi energy 
(average chemical potential of the contacts). Parame- 

.(0) 



ters of the calculation are T = 10 K, T^' = 0.005 eV 
[K = L,R),uj q = 0.05 eV, M a = 0.05 eV, j ph = 0.001 eV 
(corresponds to reorganization energy of ~ 0.05 eV). The 



source-drain voltage drop is $ = 0.01 V. Shown are self- 
consistent (solid line) and zero-order (dashed line) re- 
sults. Also shown is current profile when electron and 
phonon are decoupled (dotted line). The shift in the 
peak position is due to the reorganization energy. 




FIG. 6: Contour plot of non-equilibrium DOS vs. energy and 
position of electron level: zero-order (a) and self-consistent (b) 
results. Parameters of calculation are the same as in Fig. 

Mitra et ali^ have found that for small source-drain 
voltage ($ = ^> SOU rce - ®drain < Wo) n0 phonon side- 
band appears in the plot of source-drain current vs. 
gate volt age pot ential. This contradicts earlier findings 
(see Refs. |42|43|44|) and since consideration in Ref. 35 is 
based on the Migdal-Eliashberg theory, which is known to 
break down in the case of intermediate to strong electron- 
phonon couplingiSaLS^, this finding should be critically 
examined. The results of Fig. confirm the conclusions 
of Mitra et al. Note also that the cause of the previous 
erroneous predictions is neglect or oversimplified descrip- 



tion of hole transpor 



-59 



,42,43,44 



Finally, in the case 



where the couplings to the two contacts are proportional 
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to each other, Tl(E) = xTr(E), the current through the 
junction for resonant level model can be obtained as an 
integral over DO S 2 ^? , 

(45) 

In the case of small source-drain voltage, where T(E) is 
a smooth function of E, the current characteristics are 
determined by A(E). We can therefore demonstrate the 
"floating" behavior of phonon sidebands in a contour plot 
of the DOS vs. energy and electron level position (see 
Fig. |SJ. Parameters of the calculation here are the same 
as in Fig. EI Figure^, shows the zero-order result, while 
Figure|H|3 presents a self-consistent calculation. One sees 
that the phonon sidebands disappear around the posi- 
tion of Fermi level E — Ep — 1.95 eV. This in turn leads 
to absence of peaks in the current lineshape of Fig. [3] 
Note that since the effect is present already in the zero- 
order situation (Fig. , it can be studied analytically. 
We find, that for low temperatures (T — > 0) one can 
not see phonon sidebands in the current-voltage charac- 
teristic while changing gate voltage unless <& > ujo (see 
Appendix [EJ . 



V. CONCLUSION 

Within a non-equilibrium Green function formalism on 
the Keldysh contour, we have developed an approximate 
self-consistent procedure for treating a phonon-assisted 
resonant level model in the case of intermediate to strong 
electron-phonon interaction, where the strength of this 
interaction is determined relative to the bridge-contacts 
coupling. Our scheme goes beyond earlier considerations 
of this problem by taking into account the mutual in- 
fluence of the electron and phonon subsystems. In zero- 
order and in a single electron transport situation (Fermi 
energies in the leads far above or far below the bridge 
level so that the latter is full or empty, respectively) it 
is similar to other approachesi 2 * 4 ^^^ 4 * 4 ^. However in the 
case of a partially filled resonance level even the zero- 
order of the self-consistent scheme extends previous cal- 
culations (at least in low temperature regime) in treating 
correctly hole transport. Our approach is also similar in 
zero-order to the NLCE scheme proposed iniS, since both 
use cumulant expansion on the contour. However while 
NLCE appears to be problematic when going to the sec- 
ond cluster approximation, the present scheme is rather 
stable when higher order correlations are included by the 
self-consistent procedure. 

We have presented several numerical examples and 
compared results to those of earlier studies. The self- 
consistent calculation is found to yield drastically dif- 
ferent results as compared to the zero-order theory in 
the case of strong electron-phonon interaction and in the 
region of a partially filled electronic level. In particu- 
lar, it leads to shifts (in position and height) of peaks 



in the conductance vs. source-drain voltage plot and to 
phonon absorption signals even at low temperatures, that 
probably result from heating of the primary phonon by 
the electronic flux. The non-equilibrium electron DOS 
and the electron distribution show a similar structure, 
with peaks associated with phonon emission and absorp- 
tion. Finally, we confirm the statement of Ref. |3f| that 
current measured in a gate voltage experiment, when 
source-drain voltage is fixed at some small value, will 
not produce peaks in current vs. position of electron 
level plot, as was erroneously suggested i n 42 i 43 i 44 . Since 
the statement inS is based on application of the Migdal- 
Eliashberg theory, known to break down in the case of 
intermediate to strong electron-phonon interaction, our 
confirmation seems to be essential. Together with previ- 
ous work2l, which deals with the weak electron-phonon 
interaction case in a self-consistent manner, this provides 
tools for describing both resonant (intermediate to strong 
interaction) and off-resonant (weak interaction) tunnel- 
ing regimes in molecular junctions. 

Straightforward generalization of the present scheme 
would involve retaining mixed terms in deriving the 
equation of motion for the phonon Green function (Ap- 
pendixQ , thus abandoning the non-crossing approxima- 
tion (i.e. introducing vertex corrections). One could also 
go beyond second order in electron-phonon coupling in 
the cumulant expansion (Appendix or in coupling to 
the contacts in EOMs ( Appendices IU1 and ID)) . Develop- 
ment of a scheme spanning the entire range of parameters 
in the nonequilibrium situation is a goal for future work. 



Acknowledgments 

We are grateful to the NASA URETI program, the 
NSF-NCN program through Purdue University and the 
MolApps program of DARPA for support. The research 
of AN is supported by the Israel Science Foundation and 
by the US-Israel Binational Science Foundation. 



APPENDIX A: HAMILTONIAN 
TRANSFORMATION 

Here we derive Eq.© for the case of weak coupling 
between the primary phonon and the thermal bath (see 
below). Let us focus on the part of the Hamiltonian (|TJ 
relevant to the electron-phonon interaction 

Hel-ph = epc^c + cdoffltq + upW B bp 

+ M a Qjc + Y,U Q a Q (Al) 

13 
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We use the small polaron transformation, Eqs. @ and 
©, to get 



£o I c'c+ ujoa'a + upb^bp 

w ° ' 



(A2) 



One more transformation 

e Sa H e i- ph e- §a -> e^e^Hei-phe-^e- 8 " (A3) 

with 



(A4) 



will lead to Hamiltonian H^ t _ ph of the form (|A1|) with 
substitutions 



£ ^£o-AW and M a ^AfW (A5) 
where renormalized parameters are 

(2^) 2 



A « = M I i + y. 

W \ ^ 



/3 



(2^) 2 

"0^0 



(A6) 
(A7) 



Repeating these transformations a second time leads to a 
Hamiltoniar 
parameters 



Hamiltonian H el _ ph of the same form but with different 



A< 2 > = A< 



LO {) 



UJ LUp 



(A8) 



Wo 



(2^) 2 = 



(2C//3) 2 



LOQUJfj 



(2^/3) 



213 



E 



to 2 



2 



(A9) 



n th step 



Repeating the procedure described above, we get for the 

(2t//3) 2 



A (n) = 

i=0 



Mi") = M a 



E 



E 

(2^) 2 



UJQUJ0 



(A10) 



(All) 



Now, assuming that coupling between primary phonon 
and thermal bath is small in the sense 



E 



to 2 

OJOCO0 



< 1 



(A12) 



then continuing the procedure in the limit n — > oo leads 
to 



A (oo) = M| 



1 



Wo l-E /3 (2C// 3 ) 2 /^/3 (A13) 
= (A14) 

So we arrive at decoupling of electron and phonon degrees 
of freedom in H^} ph . Going back to the full Hamilto- 
nian 0J of the system we note that the true complete 
separation of the electronic and phononic degrees of free- 
dom, as achieved in H^} ph , is impossible here due to 
the coupling between the resonant level and the leads. 
Instead the aforementioned procedure leads to 

i7(°°) = (eo-A (oo) )c t c+ £ ^k 

+ E (^4^ (oo) +H.c.) (A15) 

ke{L,R} 

p 



where 



y(oo) v-(oo) 



exp 



M„. 



(at 



x. 



(oo) 



JJexp 

(9 



w o(l - L/3(2C / /3) 2 /wow /3 ) 

M2L^ 

^0^(1 - I] /3 (2C/ /3 ) 2 /wow / 3) 



(A16) 



(A17) 



(A18) 



Finally, neglecting in the spirit of the non-crossing ap- 
proximation the Xj°°) operators in the coupling to the 
contacts, renormalizing M a to incorporate the denom- 
inator in the exponent of X a expression, and setting 
A = A(°°), leads to ©. 



APPENDIX B: < T c X(Ti)X f (T 2 ) > IN TERMS OF 

D PaPa(T!,T 2 ) 

Here we derive Ea. l|17(l expressing the shift generator 
correlation function in terms of the phonon momentum 
Green function. We start from the Taylor series that 
expresses the correlation function as a sum of moments 



< T c X{tx)X\t 2 ) > = 



j2 f - h__)_ ( a °) r 



,1 



,1 



n—0 m—0 

X < T c Pa \n)P a m \T2) > (Bl) 
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where Eq. was used. The cumulant expansion for this 
function is 



<T c X{ Ti )X\t 2 ) >=cxp 



.p=l F ' 



= 1 +E§^i,^) 



(B2) 



P =i 



2 00 00 ^f>2 

2 E E ^[ ^Pi(n,7aV w (n,7a) 



where </? p (ti,t 2 ) is a cumulant of order p. Considering 
terms up to \ 2 a and equating same orders of A a in (|B1|) 
and (IB2I) leads to 



1 < -Pa(n) > -i < P a (T 2 ) >= Vl(n,r 2 ) 

2 < T c P (ri)P a (r 2 ) > - < P a 2 (n) > - < P a 2 (r 2 ) > 

= ^2(Tl, T 2 ) + lfil 2 (Tl,T 2 ) 

(B3) 

In steady-state < P a (n) >=< P a (t 2 ) >=< P a >, 
which implies 

Vi( t i,t 2 ) = 

Mn,r 2 ) = 2 < T c P a ( n )P a (T 2 ) > -2 < P a 2 > 

(B4) 

Using this in (|B2fl leads to (|T7jl . 

APPENDIX C: EOM FOR L> Pa P a (n, r 2 ) 

Here we derive a Dyson-like equation for the phonon 
momentum Green function Dp a p a under the Hamiltonian 
© with the EOM method. Under Hamiltonian ©, the 
EOM for the shift and momentum operators, Eqs. (0) and 
© , in the Hciscnberg picture on the Kcldysh contour are 



. dPjj) 

.dQJj) 
1 Or 



Or 
dQpjr) 



-iLU Q a {T) ~2t^2 U p Q P (t) (CI) 
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lUJ Pa(T) (C2) 

2A Q (Vkc{(T)c(T)X a (r)-R.c^ 

k£{L,R} 

~iLu p Q l3 (T)-2iU( 3 Q a (T) (C3) 

iupPp(r) (C4) 



We are looking for a Dyson-like equation for the phonon 
momentum Green function Ijl8|l . We introduce the oper- 
ator 

1 " " -~2^~o d^ +U 



Pa Pa 



(C5) 



with property 



D 



(0) -l.£)(0) 



Pa Pa 



PaPc 



(r,r')=DP aPa (r,r')-D 



^.0$) -i =S (t,t') 
(C6) 



where Dp^p is the Green function of a free phonon (de- 
coupled both from the electron and the thermal bath). 
Applying (|C5|) to H18fl from the left, taking into ac- 
count I|C1(I - I|C4() . and restricting consideration to the 
non-crossing approximation (NCA) 52 , so that terms mix- 
ing different processes are disregarded, one gets 



D { p ] p ^Dp^Ary) = S(t,t') +y / U^D PpPa (ry) 
+ i\ a J2 \Vk{~i)<T c cl(T)c{T)X a (T)P a {T r ) (C7) 

-V k *(-i) < Tj(T)Xl(T)c k (T)P a (r') > 

Next we apply operator l|C5|) to 1C7(1 from the right. The 
procedure is the same in the sense that here we deal once 
more with EOM for P a (the one depending on r'). After 
tedious but straightforward algebra and convolution of 
the result with Dp^ p we obtain a Dyson-like equation 

Dp a p a (r,r') = D^ Pa (ry) (C8) 
+ Jdn J dr 2 Df aPa (r, ri ) n PoPa (ri , r 2 ) D% Pa (r 2 , r') 



with the self-energy expression 

n PaP„( T l,T 2 ) = > — d(Tl,T 2 ) 

^ — ' Mf 



(C9) 



EKtM dp^{t 1 ,t 2 )- 1 \i y: m 



ke{L,R} 

[ 5fc (r 2 ,Ti)Gi 0) (ri,r 2 ) < T c X a { Tl )Xl{T 2 ) > +H.c. 



In what follows we neglect renormalization of the phonon 
frequency due to coupling to the thermal bath (first term 
on the right and real part of the second term in the 
self-energy expression). This is in analogy to the wide 
band approximation (for discussion on applicability of 
the approximation to coupling to the thermal bath case 
see Ref. l37|) . In the second term on the right we also 
replace lo^/loq by unity arguing that main contribution 
to Dp a p a comes from the region u>p ~ ojq. Then the 
dressed form (zero-order Green and correlation functions 
are substituted by the full ones) of (|C9|) is Eq. l(2T|) . 



APPENDIX D: EOM FOR G (ti,t 2 ) 

Here we derive a Dyson-like equation for electron 
Green function G c under Hamiltonian JfJJl using the EOM 
method. First we note that under Hamiltonian @, the 
EOA4 for operator c in the Heisenberg picture on the 
Keldysh contour is 



.dc(r) _ 
1— = £ oc(r) 



E 

ke{L,B.} 



V k *Xl(r)c k (T) (Dl) 
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with the corresponding Hermitian conjugate for the cre- 
ation operator & . We then introduce the operator 



(D2) 



with the property 
G£°> - 1 ■ ( T , T ') = (t, t') ■ G^' 1 = S(r, r') (D3) 

Here G^ (t, r') is the Green function for an electronic 
level decoupled from the contacts. Applying it to the 
electron Green function G c (t,t') first from the left and 
then from the right, one gets 

G^- 1 ■ G c (t,t') ■ G^- 1 = 5{t,t') ■ G^- 1 (D4) 
+ E \V k \ 2 g k (T,T>)<T c X a (T>)Xi(T)> 



ke{L,R} 



Finally we take convolution of l|D4|l with G^ from left 
and right and utilize Eci. (|D.3JI to arrive at a Dyson- like 
equation for G c with the dressed self-energy of the form 
presented in (25} . 



APPENDIX E: PHONON SIDEBANDS IN 
CURRENT- VOLTAGE CHARACTERISTICS 



Within the wide-band approximation, the zero-order 
lesser and greater electron Green functions (in the en- 
ergy domain) are 



G<{E) 
G>(E) 



i8(pL L - E)Y L + i6(jXR - E)T R 

{E-e y + {T/2f 
-i9(E - y, L )T L - i6{E - Li R )r R 
(E - eo) 2 + (r/2) 2 



(E3) 



(E4) 



Applying HE1(1 - HE4(I to (|T4"|) . using the resulting Green 
function in (|43|l and the spectral function A in (|45() leads 
to an expression for the source-drain current in the form 
(putting fj, L - fi R = $ > 0) 



Here we study analytically the possibility of observ- 
ing phonon sidebands in current-voltage characteristics 
when the gate voltage applied to the junction is varied. 
We restrict our consideration to the zero-order situation 
(first step of the self-consistent procedure) and low tem- 
perature (we take T = 0). The zero-order correlation 
functions for the shift generator operators (in the time 
domain) are 

< X a (t)Xl(0) > (El) 
= e -^a(2JV +i) exp { A 2[(iv + i) e -*">* + N e lUot }} 



00 \2fc 



e-^exp^e-^l^-ETT 



-ikujQt 



k=0 



< Xt(0)X a (t) > 



(E2) 



= e" A « {2No+1) exp{X 2 a [N e-^ ot + {N + l)e^ *]} 



°° y2k 



e-^explA^j-e-^E^ 



fc=0 



^_ )2k 



T - e VlVr -A 

sd ~ h r k\ 

k=0 



Mi dE 



Hr 



E_2_ 
k\ 

6{E - kuj Q - fJ, R )T R 



(E - ku - e ) 2 + (r/2) 2 

0(Hl - E - kuj )T L 
~ (E + kco - e ) 2 + (r/2) 2 



(E5) 



while the gate voltage changes the position of Eq. It is 
obvious that one can not observe phonon sidebands (k > 
terms) in the source-drain current-voltage unless $ > 

LOq. 
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